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1.  Introduction 

The  classic  machine  repair  problem  with  spares  is  a  typical 
example  of  a  finite  state  space  queueing  problem  and  consists  of  a  fixed 
number  of  identical  machines  of  which  initially  M  are  operating  and  Y 
are  spares.  The  M  machines  are  in  parallel  and  are  independent.  When 
one  fails,  it  is  instantaneously  replaced  by  a  spare  if  a  spare  is 
available;  if  not  less  than  M  machines  will  operate  until  a  repaired 
machine  becomes  aviilable.  Simultaneously,  the  failed  machine  goes 
instantaneously  into  a  repair  facility. 

1.1  Assumptions  and  Problem  Statement 

The  following  assumptions  are  made  concerning  the  machine  repair 
example. 

(a)  The  system  failure  rate  is  proportional  to  the  number  of 
operating  machines. 


(b)  Each  machine  has  exponential  failure  time  with  mean  1/A 
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(c)  There  are  c  parallel  servers  In  the  repair  facility. 

(d)  Each  server  has  exponential  service  time  with  mean  1/y  . 


Thus,  the  machine  repair  problem  is  Markovian  and  the  states  of  this 
Markov  process  can  be  described  by  a  single  number  i  ,  where  i  rep¬ 
resents  the  number  of  machines  in  the  repair  station.  The  intensity 
matrix  of  this  process  is  given  by 


Q 


~Aryi 


-x2-u2 


0 


0 

0 

0 


where 


^M+Y  -WM+Y 


MX 

> 

(0 

<  i  <  Y) 
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(Y 

<  i  <  Y+M) 

0 

9 

(i 

>  Y+M) 

i-H  . 

(0  < 

i  < 

c) 

cn  , 

(i  > 

c) 

• 

For  this  problem,  steady  state  solutions  in  closed  form  are  readily  at¬ 
tainable  [see,  for  example.  Gross  and  Harris  (1974)]. 


1.2  Transient  Solutions 

It  is  desired  to  find  the  transient  solutions  for  the  machine  re¬ 
pair  problem.  If  the  problem  has  N  states  (N  =  M+Y+l)  ,  the  inten¬ 
sity  matrix  provides  N  equations 

n'(t)  =  JI(t)  •  Q  ,  (1) 

where  H(t)  is  a  N  *  M+Y+l  component  vector  whose  elements  are 
7r^(t)  ,  the  unconditional  probability  that  the  system  is  in  state  i  at 
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> 


tine  t  ,  and  il'(t)  a  vector  of  derivatives  of  ^(t)  .  Finding  solu¬ 
tions  in  closed  fora  [see,  for  example,  Marlow  (1978)]  is  extremely  dif¬ 
ficult  and  in  most  cases  impossible  (unless  M+Y  is  very  small) . 
However,  a  variety  of  procedures  is  available  which  can  yield  numerical 
solutions  to  the  differential  equations.  In  this  paper  we  will  discuss 
several  of  these  numerical  procedures  and  attempt  to  apply  each  proce¬ 
dure  to  a  machine  repair  problem  with  the  following  parameters: 

M*4  ,  number  of  machines  initially  operating 

Y=1  ,  number  of  spares 

c=2  ,  number  of  service  channels 

A=.15  ,  machine  failure  rate 

y=.5  ,  service  rate. 


Under  the  assumptions  mentioned  in  the  preceding  section  one  can  obtain 
the  following  initialized  first  order  system  of  differential  equations. 
Assuming  at  t=0  all  machines  are  up: 
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with  initial  value  11(0)  =  [1,0, 0,0, 0,0]  . 


2.  Randomization  Techniques 

The  randomization  procedures  give  a  method  of  calculating  the 
transition  probability  matrix  P(t)  ,  i.e.,  the  order  (M+Y+l)  square 
matrix  whose  elements  are  p^  (t)  ,  the  probability  that  the  system  is 

in  state  j  at  time  t  given  that  it  started  in  state  i  (0  <  i,j  < 
M+Y)  .  From  P(t)  ,  by  using  the  initial  probability  vector  11(0)  , 


-  3  - 


T-436 


n(t)  can  be  calculated  from  IT ( t)  -  Il(0)P(t)  .  We  denote  the  (i,j )th 
element  of  Q  by  q^  ,  and  define  a  scalar  3  as 


8  *  max  |q 


ii1 


(2) 


Now  define  a  matrix  P  with  elements  as 

p  *  {pij}  5  i  +  ih  • 


where  1  is  an  identity  matrix  having  the  same  order  as  Q  . 


The  matrix  P  is  stochastic,  and  it  has  been  shown  [see  Cohen 
(1969)]  that  the  elements  of  P(t)  can  be  obtained  by 


00 


I 

n=0 


(3) 


where  *s  the  (i,j)£?z  element  of  matrix  P  raised  to  the  nth 

power.  This  method  is  called  randomization  because  it  can  be  inter¬ 
preted  as  a  discrete  time  Markov  chain  with  transition  probabilities 
p^  and  transition  time  generated  by  a  Poisson  process  at  rate  8  . 


To  compute  p^ (t)  involves  raising  the  matrix  P  to  the  nth 

power  for  n=0,l,2,...  .  The  numerical  procedure  truncates  n  at  some 
appropriate  value,  say  m  ,  so 


p. .  (t)  = 


-3t 


n-l 

l 

i-O 


8ntn 


.(n) 


n!  rij 


+  R 


m 


(4) 


where  R  is  the  error  due  to  truncation, 
m 


2.1  Barzily  and  Gross  Method 

Barzily  and  Gross  (1979)  proposed  a  criterion  to  truncate  n 

such  that  for  a  given  e>0  ,  the  smallest  m  is  chosen  such  that  the 

error  R  obeys 
m 
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R  < 

m  = 


1 


e-3 K  e 


n=0 


n! 


(5) 


Their  paper  shows  that  such  an  m  exists.  For  a  machine  repair  problem 

_3 

with  parameters  M  ■  Y  =  C  =  1  and  £  =  10  ,  the  result  of  applying 

this  criterion  has  been  shown  to  be  quite  satisfactory,  in  fact  accurate 
up  to  the  second  decimal  place.  A  complete  description  of  the  Barzily- 
Gross  method  together  with  the  study  of  the  transient  effects  and  the 
speed  of  convergence  to  steady  state  for  machine  repair  problems  can  be 
found  in  the  above  cited  reference.  The  algorithm  has  been  coded  in 
FORTRAN  to  run  on  the  IBM  3031  computer,  under  the  program  name  WONG. 

The  following  section  is  a  brief  description  of  program  WONG. 


2 . 2  Program  WONG 

WONG,  a  FORTRAN  code  originally  designed  to  find  Lht  spares  in¬ 
ventory  level  and  number  of  repair  channels  necessary  to  guarantee  a 
prespecified  service  level  for  a  machine  repair  problem,  included  in  it 
a  program  to  provide  transient  solutions  for  the  system  state  probabil¬ 
ities.  This  portion  of  the  program  was  separated  from  the  original  and 
updated  to  stand  alone  as  a  provider  of  transient  solutions  to  machine 
repair  problems.  Input  requirements  for  program  WONG  are  given  in 
Appendix  1.  Results  of  applying  program  WONG  to  the  sample  problem  of 
Section  1.2  are  given  in  Section  5. 


3.  The  QUE  Package 

Grassmann  and  Servranckx  (1979)  developed  a  FORTRAN  based  package 
for  finding  transient  solutions  for  moderate  sized  queueing  networks  (up 
to  ten  state  variables).  The  method  adopted  in  this  package  is  in  fact 
based  on  the  randomization  procedure  discussed  in  the  preceding  sections. 
The  truncation  criteria  are  somewhat  different  and  are  described  in  the 
following  section.  We  have  adopted  the  sample  problem  of  Section  1.2  to 
the  specifications  of  this  package  and  the  results  are  presented,  along 
with  those  of  program  WONG,  in  Section  5. 
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Grassmann  (1977)  has  also  shown  that  the  truncation  error,  R 


m 


can  be  bounded  with  reasonable  accuracy;  that  Is, 

R  <  1 


m 


m-1 

I 

n=0 


(et)n  e-gt 


nl 


(6) 


For  small  0t  ,  the  sum  in  the  RHS  of  (6)  can  easilv  be  evaluated  for 

fixed  m  ,  and  thus  m  can  be  determined  such  that  R  will  be  below 

m 

a  prescribed  value  e  >  0  ,  as  was  done  by  Barzily  and  Gross  (1979).  For 
large  0t  one  can  approximate  the  Poisson  distribution  by  the  normal 
distribution.  In  the  QUE  package,  m  is  set  equal  to 

m  ■*  0t  +  4/0t  +  5  .  (7) 

Using  Poisson  tables  for  0t  £  20  ,  or  normal  tables  otherwise  this 

-  4 

procedure  guarantees  that  such  an  m  yields  an  R^  less  than  10 


3.1  Formulation  of  the  Sample  Problem 
for  the  QUE  Package 

To  solve  the  machine  repair  problem  by  using  the  QUE  package, 
one  must  formulate  the  problem  to  fit  the  network  structure  input 
requirements.  The  output  statistics  are  then  obtained  from  the  program. 
One  way  of  formulating  the  problem  to  fit  the  QUE  package  requirements  is 
shown  in  Figure  1. 


Operating  units 

Repair 

Spares  inventory 

r 

faci 

lity 

f 

Figure  1. — Formulation  of  the  sample 
machine  repair  problem  for 
the  QUE  package. 
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It  is  necessary  to  define  state  variables  and  event  descriptions 
together  with  type,  conditions,  net  effect,  and  rate  parameters.  These 
are  shown  in  Table  1. 


TABLE  1 

THE  EVENT  DESCRIPTION  OF  THE  PROBLEM  FOR  THE  QUE  PROGRAM 


Events 

Type 

Rate 

Condition 

Net  Effect 

1.  Arrival  into  repair 
station 

a.  when  no  machine 
is  in  repair 
station 

1 

.6 

X1=0 

(+1) 

b.  when  the’ e  are 
some  machines 
in  the  repair 
station 

1 

.75  -  .15X1 

1  <  Xi  <  4 

(+1) 

2.  Service:  when  only 
one  repairman  is 
busy 

1 

.5 

Xl-1 

(-1) 

3.  Service:  when  both 
repairmen  are  busy 

1 

1 

2  <  X!  <  5 

(-1) 

The  state  variable  is  described  by  XI  =  number  of  machines  in 
repair  station,  0  <  XI  <_  5  . 

Each  event  has  a  rate  function  which  associates  the  rate  of  each 
transition  with  the  starting  state.  The  rate  functior  may  be  constant  or 
a  function  of  the  state  variable.  Types  of  events  are  classified  as  type 
one,  having  finite  rate  event,  or  as  type  two,  having  infinite  rate 
event.  The  state  space  of  the  system  is  defined  by  general  conditions 
represented  by  linear  inequalities  involving  the  state  variables.  The 
net  effect  is  the  \alue  of  the  state  variable  after  the  event  occurs  and 
is  determined  by  incrementing  or  decrementing  the  value  by  a  constant 
prior  to  the  event’s  occurrence.  A  more  detailed  explanation  can  be 
found  in  Grassmann  and  Servanckx  (1979).  Appendix  2  shows  the  QUE  pro¬ 
gram  input  requirements  for  the  sample  problem. 
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4.  Numerical  Integration  Methods 


Numerical  integration  methods  can  be  employed  to  solve  a  system 
of  ordinary  differential  equations  described  by 


”y{(t)" 

f1(y1,...,yp;t) 

y'2U) 

II 

f2(y1,...,yp;t) 

fp(Yi»  * • • »yp;t) 

f(Y,t) 


(8) 


with  known  initial  value  Y(tg)  •  The  standard  techniques  are  generally 

variations  of  either  Runge-Kutta  (R-K)  or  predictor-corrector  (p-c)  methods. 

Runge-Kutta  methods  are  based  on  formulas  that  approximate  the 
Taylor  series  solutions 


2  k 

yi(t-Hi)  -  yi(t)  +  h*y|(t)  +  ~  y'^(t)  +  ...  +  £7  y^(t>  +  » 

i=l,...p  .  These  methods  use  approximations  for  the  "second  and  higher- 
order  derivatives.  Euler's  method  is  a  special  R-K  method,  with  Jt=l  . 
These  methods  have  been  used  by  several  authors  [e.g.,  Bookbinder  and 
Kartell  (1379) ,  Grissmann  (1977),  Liitschwager  and  Ames  (1975)  and  Neuts 
(1975)]  to  find  transient  solutions  in  queueing  systems. 


The  predictor-corrector  methods  require  information  about  several 
previous  points  in  order  to  evaluate  the  next  point.  These  methods  in¬ 
volve  using  one  formula  to  predict  the  next  Y(t)  value,  followed  by  the 
application  of  a  more  accurate  corrector  formula.  Unlike  the  R-K 
methods,  p-c  methods  are  not  self-starting;  hence,  they  use  the  R-K 
method  to  obtain  the  first  Y(t)  value. 


Predictor-corrector  methods  can  provide  an  estimate  of  the  local 
truncation  error  at  each  step  in  the  calculations,  in  contrast  to  the  R-K 
methods,  which  cannot  obtain  such  an  estimate. 
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Predicto c-corrector  methods  include  Milne’s  method,  Hamming's 
method,  and  Adams'  methods.  Methods  based  on  the  Adams  formulas  have 
performed  very  well  in  test  problems  [see  Hull,  Enright,  Fellen,  and 

Sedgwick  (1972)]  even  for  nonstiff ^  systems  or  when  function  evaluations 
are  relatively  expensive.  Hull  et  al.  also  concluded  that  R-K  methods 
are  not  competitive,  although  fourth  or  fifth  order  methods  are  best  for 
problems  in  which  function  evaluations  are  not  very  expensive  and  accura¬ 
cy  requirements  are  not  very  stringent. 

Predictor-  corrector  methods  have  been  used  by  Ashour  and  Jha  (1973) 
for  queueing  problems. 

A  variety  of  routines  is  available  for  solving  a  system  of  ordi¬ 
nary  differential  equations.  They  include  RKGS,  DRKGS  (fourth  order  R-K 
formulas),  HPCG,  DHPCG,  HPCL,  DHPCL  (Hamming's  Methods),  all  from  the  IBM 
Scientific  Subroutine  Package  [IBM  (1968)];  and  DVERK  (Verner's  fifth  and 
sixth  order  R-K  formulas)  and  DVOGER  (Gear's  Method)  in  the  International 
Mathematical  and  Statistical  Libraries  (IMSL)  package  [IMSL  (Ed.  6)].  One 
routine  based  on  extrapolation  methods  is  DREBS,  also  in  the  IMSL 
package,  which  uses  the  Bulvisch-Stoer  method. 

4.1  Gear's  Algorithm 

C.  W.  Gear  (1971a,  1971b)  proposed  a  variable-order  integration 
method  based  on  Adams'  predictor-corrector  formulas  of  orders  one  through 
seven.  It  uses  an  order  one  formula  to  start  and,  for  this  reason,  must 
start  with  very  small  step  size  when  the  error  tolerance  is  stringent. 

Gear's  algorithm  includes  a  special  approach  for  dealing  with 
stiff  differential  equations. 


A  stiff  system  of  ordinary  differential  equations  is  character¬ 
ized  by  the  property  that  the  ratio  of  the  largest  to  the  smallest 
eigenvalue  is  much  greater  than  one. 
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The  Adame'  formulas  fall  under  two  general  categories — open  and 
closed  formulas .  The  Adams-Bashforth  pth  order  formulas  (open  formula) 
can  be  written  as 


Vl  +  h  > 


(9) 


where  y±  =  y(t1>  ,  t± 


ih  ,  y^  =  f(y^,t^)  .  The  order  of  the  method 

is  one  less  than  the  order  of  the  truncation  error  per  step.  The  Adams- 
Moulton  pth  order  formulas  (closed  formulas)  can  be  written  as 


P-1 


n,m 


»X-K 


(10) 


The  coefficients  3^  and  6*  are  given  by  Henrici  (1962) .  Equation  (9) 

is  used  as  the  first  approximation  in  Equation  (10).  Thus  (9)  is  used  as 
the  predictor  eqt ation  and  (10)  as  the  corrector  equation.  Whenever  (10) 
converges  (as  is  true  when  h  is  small  and  f  is  smooth)  ,  the  trunca¬ 
tion  error  introduced  at  the  nth  integration  step  is 

h1*4^-  y^P4"!)  (t^)  +  0(h**^)  ,  where  y^  is  the  kth  derivative  of 

y  ,  and  are  constants  [see  Henrici  (1962)]. 


The  predictor  equation  (9)  is  equivalent  to  fitting  a  pth  degree 

polynomial  through  the  known  quantitites  y  .  ,  hy'  ,,  ...»  hy' 

n-i  n-I  n-p 

For  more  details  of  the  algorithm,  see  Gear  (1971a,  1971b). 


4.2  DVOGER  Subroutine 

DVOGER  is  a  FORTRAN  routine  based  on  Gear's  algorithm  designed  to 
solve  a  set  of  first  order  ordinary  differential  equations.  The  algo¬ 
rithm  chooses  the  order  of  approximation  such  that  the  step  size  is  in¬ 
creased,  thereby  decreasing  the  solution  time.  The  option  of  using  a 
particular  method  is  done  through  a  switch  variable  (MTH) .  Results  of 
using  DVOGER  on  the  sample  problem  are  also  given  in  section  5.  Appendix 
3  shows  the  input  and  programming  requirements  for  exercising  DVOGER. 
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Conclusions  and  Numerical  Results 


We  have  nresented  three  numerical  methods  in  computing  the  tran¬ 
sient  probabilities  for  finite  Markovian  queues.  Transient  probabilities 
often  provide  a  realistic  picture  of  actual  queueing  systems,  and  some¬ 
times  it  is  desirable  to  know  how  fast  they  converge  to  steady  state 
[Barzily  and  Gross  (1979)]. 

The  methods  considered  in  this  paper  fall  into  two  categories, 
namely,  randomization  and  predictor-corrector  numerical  integration.  In 
general  these  methods  give  reasonably  accurate  results  for  a  moderate 
sized  problem.  Table  2  shows  the  output  of  these  programs  for 
t*l»3,5,7,9,12.  The  QUE  package  and  DVOGER  show  almost  equal  results, 

while  program  WONG  deviates  from  the  other  two  by  at  most  3x10 ,  which 
is  reasonably  compatible. 

In  terms  of  set-up  effort,  program  WONG  gives  the  least  degree  of 
difficulty  since  it  was  written  primarily  for  machine  repair  problem. 

The  biggest  concern  with  respect  to  the  QUE  program  was  the  huge  core 
storage  requirement,  which  exceeds  the  current  daytime  capacity  of  384K 
bytes  of  the  IBM  3031  at  the  GWU  Center  for  Academic  and  Administrative 
Computing.  Future  modification  by  redimensioning  is  suggested.  In  using 
DVOGER,  one  must  carefully  choose  applicable  parameter  values,  as  in  the 
step  size.  Total  running  times  of  the  programs  are  2.35  seconds  for  pro¬ 
gram  QUE,  6.13  seconds  for  program  WONG,  and  162.62  seconds  for  DVOGER. 
The  reason  for  the  length  of  the  latter  is  that  with  step  size  fixed  at 
-4 

1x10  ,  DVOGER  mist  be  called  120,000  times  to  integrate  for  each  time 

point  from  t*»0  to  t*12  .  There  is  a  need  to  explore  further  the  best 
options  of  DVOGER  to  find  those  which  might  reduce  running  time  consider- 
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TABLE  2 

PROGRAM  OUTPUT  OF  WONG,  QUE,  DVOGER  FOR  THE  TRANSIENT  SOLUTION 
OF  A  MACHINE  REPAIR  PROBLEM  WITH  PARAMETERS  M-4,  Y-I,  C-2, 


X“.15,  p*.5 

and  11(0) 

-d.o.o 

,0,0,0) 

AT  t-1,3 

,5,7,9, 

12 

Time 

Program 

Vt} 

7^(0 

7T2(t) 

7r3(t) 

Vc> 

7t5(t> 

1 

WONG 

.6235 

.2946 

.0708 

.0096 

.0007 

.0000 

QUE 

.6237 

.2949 

.0709 

.0097 

.0007 

.0000 

DVOGER 

.6237 

.2949 

.0709 

.0097 

.0007 

.0000 

3 

WONG 

.3944 

.3686 

.1722 

.0536 

.0099 

.0008 

QUE 

.39  4  5 

.3688 

.1723 

.0536 

.0099 

.0008 

DVOGER 

.3945 

.3688 

.1723 

.0536 

.0099 

.0008 

5 

WONG 

.3331 

.3663 

.2005 

.0781 

.0192 

.0022 

QUE 

.3333 

.3665 

.2006 

.0782 

.0192 

.0022 

DVOGER 

.3333 

.3665 

.2006 

.0782 

.0192 

.0022 

7 

WONG 

.3119 

.3618 

.2091 

.0886 

.0244 

.0033 

QUE 

.3122 

.3620 

.2093 

.0887 

.0244 

.0033 

DVOGER 

.3122 

.3621 

.2093 

.0887 

.0244 

.0033 

9 

WONG 

.3038 

.3595 

.2123 

.0931 

.0269 

.0038 

QUE 

.3039 

.3597 

.2124 

.0932 

.0269 

.0038 

DVOGER 

.3039 

.3597 

.2124 

.0932 

.0269 

.0038 

12 

WONG 

.2995 

.3580 

.2138 

.0955 

.0283 

.0042 

QUE 

.2997 

.3582 

.2140 

.0956 

.0284 

.0042 

DVOGER 

.2997 

.3582 

.2140 

.0956 

.0284 

.0042 
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APPENDIX  1:  PROGRAM  WONG 


A.  Input  Requirements 


TABLE  Al.l 

PARAMETER  CARD  INPUT 


Columns 

Format 

Input  Name 

Explanation 

1-5 

15 

M 

Number  of  machines  initially  operating 

6-10 

15 

IC 

Number  of  service  channels 

11-15 

15 

IY 

Number  of  spares 

16-27 

F12.7 

RLAM 

Machine  failure  rate  (Poisson  mean) 

28-39 

F12./ 

RMU 

Service  rate  (Poisson  mean) 

40-49 

F10.6 

EPS 

Tolerance  value 

TABLE  AI.2 
TIME  INPUT 


Columns  Format 

Input  Name 

Explanation* 

mm 

TDEL 

Time  T 

*Time  at  which  transient  probability  is 
required  is  format  free,  but  it  must  be  coded 
starting  from  column  one,  and  a  separate  C3rd 
is  required  for  every  time  desired. 

B.  Numerical  Example  Input 

We  shall  illustrate  the  use  of  program  WONG  on  the  sample  problem 

-3 

given  in  Section  1.2.  The  cards  for  the  sample  problem,  with  £  *  10 
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36 


and  t  =  3,5  are  shown  in  FigureAl.l.  The  output  obtained  for  this  prob¬ 
lem  is  tabulated  in  Table  A1.3. 


Card  Type 

Card  Image 

Job 

//  STANDARD  JOB  CARD 

JCL 

//  EXEC  F0RT2 

JCL 

//F0RT.SYSIN  DD  * 

Program  Deck 

[Program  WONG  deck) 

JCL 

//Gtf.SYSIN  DD  * 

Parameter 

A  2  1  .15  .5  .001 

Time 

3 

Time 

5 

JCL 

// 

Figuie  Al.l — Card  input  program  WONG  foi  sample 
problem. 


TABLE  A1.3 


THE  OUTPUT  OF  PROGRAM  WONG 
FOR  THE  SAMPLE  PROBLEM* 


■ 

n(t) 

■ 

=» 
o 
.  ^ 

rr 

^(t) 

iT2(t)  ir3(t) 

Ln 

/•s 

r* 

3 

.3686 

.1722  .0536 

.0099 

.0008 

5 

.3663 

.2005  .0781 

.0192 

.0022 

*Note 

:  the 

initial 

distribution  of 

the  system 

is 

assumed  to  be  tt0(0)-1,  tt^O-O,  i>l  . 
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APPENDIX  2:  PROGRAM  QUE 


Input  Deck  Format  Specification 


TABLE  A2.1 


INPUT  DECK  FORMAT  FOR  PROGRAM  QUE 


Problem  title  card: 

Function:  for  documentation  only 


Columns  Format 


Field  description 


Problem  title 


2.  Problem  specification  card: 

Function:  describes  the  number  of  state  variables,  number  of 
events,  and  the  number  of  general  system  conditions 


Columns  Format 


Field  description 


L-2  12  Number  of  state  variables 

)-4  12  Number  of  events 

i-6  12  Number  of  state  space  restrictions 

3.  Maximum  vector  card  (one  card  for  each  state  variable — for 
machine  repair  problems  only  one  is  required) 

Function:  to  describe  the  highest  possible  value  of  each  state 
variable  (a  maximum  of  ten  state  variables) 


Columns  Format 


Field  description 


L-2  12  Maximum  value  of  state  variable  1 

4.  Event  title  cards  (one  for  each  event) 

Function:  for  documentation  purposes  only 


Columns  Format 


Field  description 


20A4  Event  title 
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TABLE  A2.1 — continued 


5.  Event  specification  card  (one  for  each  event) 

Function:  indicates  type  of  event  and  the  rate  of  this  event 


Columns 

Forma* 

Field  description 

1 

11 

Print  flag  for  transitions  (1  «*  Yes) 

2 

11 

Event  type  (1  for  machine  repair  problem) 

3 

11 

Number  of  specific  conditions  (0  *or  machine  repair 
prob lem) 

4 

11 

Flag  for  new  minima  and  maxima  of  the  state  variable 
(1  -  Yes) 

5-10 

F6.2 

Rate  of  this  event 

11-12 

12 

State  variable  on  which  rate  depends  (if  zero,  rate 
is  a  constant) 

13-18 

F6.2 

Increase  of  rate 

6.  New  maxima  and  minima  vector  card  (one  for  each  event) 

Function:  resets  the  maximum  and  minimum  values  for  the  state 
variable 

Columns  Format _ Field  description _ 

1-2  12  New  maximum  for  state  variable  XI 

3-4  12  New  minimum  for  state  variable  XI 

7.  Net  effect  ''ard 

Function:  defines  the  function  f(x)  which  converts  the  starting 
state  into  the  target  state 

Columns  Format _ Field  description _ 

1-3  13  Net  effect  for  state  variable  XI 

8.  Trailer  card 

Function:  delimiter  card  used  to  indicate  the  end  of  events 
section  of  the  system's  input.  The  first  four  bytes  of  the 
record  must  contain  the  string  "END" 

Columns  Format _ Field  description _ 

1-4  A4  Control  field  (value  is  "END") 

5-80  19 A4  Ignored 


TABLE  A2.1 — continued 


9.  Probability  specification  card 

Function:  gives  the  number  of  nonzero  initial  probabilities 
(only  the  nonzero  ones  need  to  be  entered) 


Columns 

Format 

Field  description 

1-2 

12 

Number  of  nonzero  initial  probabilities 

3-8 

F6.2 

The  starting  time  of  the  system 

10. 

Initial  probability  card 

Function:  specified  initial  probability  and  the  state  to  which 
it  pertains  (one  for  each  state  variable — only  one  required  for 
machine  repair  problem) 

Columns 

Format 

Field  description 

1-6 

F6.5 

Initial  probability 

7-8 

12 

State  variable  XI 

11. 

Time  specification  card 

Function:  to  indicate  the  time  for  which  transient  solutions  are 
required,  to  indicate  what  measures  are  to  be  printed,  and  to 
give  a  criterion  whether  or  not  to  continue  calculating  the 
times  on  the  following  card 

Columns 

Format 

Field  description 

1 

11 

Number  of  times  for  which  solutions  are  required 
(5  or  less) 

2 

11 

Print  flag;  if  value  is  1,  joint  distributions  are 
printed 

3 

11 

Print  flag;  if  value  is  1,  marginal  distributions 
are  printed 

4 

11 

Print  flag;  if  value  is  1,  expectations  are  printed 

5 

11 

If  value  is  1,  other  cards  follow  (having  the  same 
format  as  this  one) 

6-11 

F6 . 5 

Stopping  criterion  (accuracy  desired) 

12. 

Time  card 

Function:  gives  each  time  (a  maximum  of  five)  for  which  results 
are  desired 


A 


T-436 


TABLE  A2.1 — continued 


Columns  Format 


1-5  F5.2  First  time  t. 


6-10  F5.2  Second  time  t„ 


Field  description 


the  times  t^,t_,...  must  be  input  in  increasing  order  of  magnitude 


B.  Numerical  Exanple  Input 


A2.2. 


Input  data  for  the  sample  machine  repair  problem  is  shown  in  Table 


TABLE  A2.2 


INPUT  CARDS  FOR  PROGRAM  QUE 


Card  Input 


TRANSIENT  SOLUTION  FOR  MACHINE  REPAIR  PROBLEM 

010400 

05 

ARRIVAL  INTO  REPAIR  STN:  (A)  NO  MACHINE 
1101000.60 


ARRIVAL  INTO  REPAIR  STN:  (B)  1  MACHINE 
1101000.7501-00.15 


SERVICE  WHEN  ONE  REPAIRMAN  IS  IDLE 
1101000.50 
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TABLE  A2.2 — continued. 


|  — i  ■■  r- 

Card 
Type 


Card  Input 


4  SERVICE  WHEN  BOTH  REPAIRMEN  ARE  BUSY 

5  1101001.00 

6  0502 

7  -01 

8  END  OF  EVENT  SECTION 

9  01000.00 

10  1.000000 

11  41111.0001 

12  00.2500.5000.7501.00 

11  41111.0001 

12  01.2501.5001.7502.00 

11  41111.0001 

12  02.2502.5002.7503.00 

11  41111.0001 

12  03.2503.5003.7504.00 

11  *#1111. 0001 

12  04.2504.5004.7505.00 

11  41111.0001 

12  05.2505.5005.7506.00 

11  41111.0001 

12  06.2506.5006.7507.00 

11  41111.0001 

12  07.2507.5007.7508.00 

11  41111.0001 

12  08.2508.5008.7509.00 

11  41111.0001 

12  09.2509.5009.7510.00 

11  41111.0001 

12  10.2510.5010.7511.00 

11  41 111.0001 

12  11.2511.5011.7512.00 
11  00000 
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APPENDIX  3:  DVOGER  SUBROUTINE 


A.  Input  and  Options 

This  section  discusses  the  input  requirements  for  DVOGER.  Since 
DVOGER  is  a  library  subroutine,  one  must  write  a  computer  program  in 
order  to  use  it.  The  advantage  is  in  the  flexibility  of  inputting  the 
parameter  values,  as  well  as  in  the  choice  of  output  variables,  frequency 
of  printing  the  solutions,  and  so  forth. 

The  input  structure  is  as  follows. 

1.  Job  and  JCL  cards.  See  Section  B  for  the  standard  and  job 
control  language  cards. 

2.  Main  program.  The  main  or  the  calling  program  is  to  be 
written  in  FORTRAN.  The  proper  dimensioning  of  arrays,  the 
input  mode  of  parameter  values,  the  nunfeer  of  calls  to  DVOGER, 
and  the  frequency  of  printing  the  solution  must  be  determined 
by  the  user.  Moreover,  an  external  subroutine  DFUN  is  to  be 
written  by  the  user  to  compute  functional  values  F(y,t)  or 
the  Jacobian  of  F(y,t)  . 

The  parameters  needed  for  the  main  program  include: 

N  *  number  of  first  order  differential  equations 

M  *  order  of  Jacobian  (M«N) 

T  *  initial  value  of  independent  variable  (e.g.,  time) 
MTH  *  method  indicator 

0,  predictor-corrector  (Adams)  method 

1,  variable-order  method,  suitable  for  stiff 

-  systems  (partial  derivatives  provided  by  user) 

2,  variable-order  method  (partial  derivatives  com¬ 
puted  by  numerical  differencing) 
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Y(1,N) 

YMAX(N) 

HMIN 

HMAX 

H 


JSTARI 

EPS 


an  input  array  of  initial  solutions  at  T  ;  array 
Y  is  to  be  declared  8xN 

an  input  array  of  maximum  absolute  value  of  solu¬ 
tion 

smallest  step  size  allowable 
maximum  step  size  allowable 

step  size  to  be  attempted  on  the  next  step;  this 
is  to  be  used  if  it  does  not  cause  a  larger  error 
than  requested 

-1,  repeat  the  last  step  with  a  new  H 

0,  initialize  the  integration  (for  first  call  to 
DVOGER) 

1,  take  a  new  step  continuing  from  the  last 

maximum  error  criterion  such  that  the  single  step 
error  estimates  divided  by  YMAX(l)  are  less  than 
EPS  in  norm. 


The  call  to  DVOGER  is  done  by  the  statement, 

CALL  Dv0GER  (DFUN,  Y,  T,  N,  MTH,  MAXDER,  JSTART,  H,  HMIN, 
HMAX,  EPS,  YMAX,  ERROR,  WK,  IER) . 


3.  Subroutine  DFUN.  DFUN  is  user-supplied  and  is  to  be  declared 
by  an  EXTERNAL  DFUN  statement  in  the  main  or  calling  program. 
DFUN  specifies  the  problem  for  DVOGER.  It  provides  the  system 
of  equations  and  the  Jacobian.  The  parameters  Include: 


YP(1,N>  ■  vector  of  solution  TP 

TP  *  present  time 
M  *  order  of  Jacobian 

0,  DFUN  computes  F(YP,TP) 

1,  DFUN  computes  Jacobian  of  F  evaluated  at 
(YP.TP) 


YP  is  to  be  declared  as  an  8*n  array. 


Numerical  Exaim  le  Input 


We  shall  illustrate  the  solution  of  the  sample  problem  using  the 
DVOGER  subroutine.  The  structure  of  a  single  Job  run  is  as  follows. 
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(a)  Job  and  JCL  cards 

//  STANDARD  JOB  CARD 
//  EXEC  F0RC2 

//F0RT.SYSIN  DD  * 

[Main  Program  and  Subroutine] 

//G0.SYSLIB  DD 

//  DD  DSN=GWU.IMSL6 .LM0D .D,DISP=SHR 

//  DD  DSN“GWU . IMSL6 . LM0D . S , DISP*SHR 

//G0.SYSIN  DD  * 

(b)  Figure  A3. 1  shows  the  main  program  and  subroutine  DFUN,  and  the 
input  and  our  option s,  where  the  transient  solutions  are  re¬ 
quired  for  times  T=l,2,  and  3  . 

(c)  Program  Output.  The  step  size  h  was  fixed  at  0.0001  by 
specifying  HMAX^HMIN™!!^  .0001  .  The  method  used  was  the 
predictor-corrector  method  based  on  the  Adams  formulas 
(MTH-0)  .  The  output  of  transient  solutions  was  printed  out 
at  every  time  increment  of  .01  starting  at  t=0  to  t=3  . 

-7 

The  error  tolerance  was  set  at  10 

The  program  below  was  written  for  this  specific  problem,  although 
a  more  general  program  can  be  written  to  handle  any  problem  with  arbi¬ 
trary  parameter  values. 

The  authors  have  not  tested  this  subroutine  for  options  which  give 
minimum  execution  time,  as  this  was  not  the  purpose  here. 
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DIMENSION  YP!6, 10) ,PW(  10.  I0).Y!8,  1 0) , YMAX ( 1 0) ,DY!  10) 

DIMENSION  NK! 140,20), ERROR! I0),C! 10) 

UOUDLc  PRECISION  YP 

DOUBLE  PRECISION  C,  Q,  E,H, P,  T, Y, Rl  ,  BND, EPS ,  EUP,EI)WM,  EM0 1  ,01 .0’, 

*  EN02 , ENQ3, HM AX , HM I M , HN Erf .HOLD, TOLD , YMAX , 

*  ERROR,  RACUM.WK.XK,  ZERO,  HALF.  ONE,  OMRP.P.U 
EXTERNAL  DF'jN 

KK=0 

N«6 

M=6 

T=0.0i)0 
Y(  1 , 6 ) =0 . ODO 
Y  ( I  , 5 ) =0 . 000 
Y ( I , 4 ) =0 . 000 
Y< l , 3)«0.0D0 
Y( I ,2)=O.ODO 
Y( I , I )=l .000 
YMAX ( 1 )=J .000 
YMAX(2)*J  .OEO 
YMAXC3)=I .000 
YMAX ( 4 )= I .000 
YMAX(5)=I .000 
YMAX(6)=1 .000 
JSTA RT=0 
IN0=0 
MTri=0 

HMAX=I .OD-4 
HMIM=I .00-4 
'  H* 1 .00-4 
EPS»l.0D-7 
WRITE! 6, 1  00 ) 

100  FORMAT! '0',9X, 'T' ♦ l 2X, 'PO' , J2X,/PIy, 12X,'P2', I2X,'P3', 

*  12X,/P4-/,  12X,  /p'5*  ) 

WRITE! 6,200)  T,  Y!  1 ,  I  ) ,  Y<  I  ,2).  Y(  1 , 3),  Y(  I  .  4  ) .  Y(  I  ,5  )  ,Y(  1  .6) 

DO  JO  K=1 ,30000 

CALL  DV0GER!DFUN, Y,T,N, WTH.MAXDER, JSTART.H.H1IN. HMAX, EPS, 

*  YMAX,  ERROR,  WK, IER) 

KK=  KK+ 1  ■ 

If IKK. NE. 100)  GOTO  250 

WRITE! 6,200)  T, Y! 1 , 1 ) , Y! 1 ,2) , Y! 1 , 3) , Y! I ,4 ) , Y! 1 ,5 ) , Y! 1 , 6 ) 

200  FORMAT! '0',7!2X,F10.8)) 

KK=*0 

250  CONTINUE 
H=1 .00-2 
IND-0 
MTH=0 

IF(MTH.NE.I)  GOTO  10 
PW( I , l >=-.6 
PW( 1,2)-. 5 
PWC I ,3 )«0.0 

Figure  A3. 1— Program  listing  to  call  DVOGER  for  the  machine  repair 
problem,  t-1,2  and  3  . 
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Ph<  I  ,4)=Q,0 
Pn( I , 5 )=0.0 
PW(  I  ,6)=0.G 
PW(2,  I  )  =  .  6 
PW(2,2>=-1.  1 
PN( 2,3 )= I  .0 
Ph ( 2 , 4 ) =0.0 
PN<2,5)=0.0 
P<«  ( 2 , 6 )  =0 .0 
PW( 3, I )=0.0 
P(l(3,2)  =  .6 
Prt(3,3)*-I.45 
P 3 , 4 )  =  I  .0 
Prt(3,5)=0.0 
Pw(3,6)=0.0  . 

P  rt  ( 4 , I )=0.0 
P<>  (4,2  )=0,0 
P-.C  4,3)=. 45 
P.<(  4,4  )=-  | .  3 
P><(  4,5)=l  .0 
P»i  (  4, 6)  =0,0 
PW( 5, i >=G„0 
Prt ( 5, 2 )=0.C 
PH (5,3 )=0.0 
PN(5,4)=.3 
PW ( 5, 5 )=-  I . >5 
PW  (  5, 6  )  =  I  .0 
P/K6,  I  >*0.0 
PH ( 6, 2 )=0. 0 
PW ( 6,3 ) =0,0 
PH ( 6,4 ) =0,0 
PN<6,5)=. 15 
PW( 6,6 )=- 1 ,0 
io  continue 

STOP 

END 

SUBROUTINE  DFUNI YP,TP,M,DY,PW, IND) 

DIMENSION  PHI  1 0, 1 0) , YP( 8, 1 0) , DY( 1 0) 

DOUBLE  PRECISION  YP,TP, DY, PH  * 

IF(IND.EO.O)  GOTO  5 

PW(  1,1 )=-,6 

PW( I  ,2)=,5 

PW( 1 ,3)=0.0 

PW( I ,4  >=0.0 

PH(  I , 5)=0.0 

Pn(  1 ,6)=0,0 

P W(2,  I  )  =  .  6 

PN(2,2 )-- I . I 

Pw( 2, 3 )=l .0 

PW(2,4)=0.0 

PW<2,5>=0.0 

PN (2,6 )=0.0 

PH (3, I )=0.0 

PW ( 3 , 2 ) * . 6 

Prt(3,3)=-I,45 

Figure  A3 . 1 — continued 
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PW<3,4)=I  .( 

Prt(  3 ,5 ) -*0.0 
PW( 3,6 )=0,0 
PW (4, ] )=0.0 
PW( 4,2 )=0.0 
Prt( 4,3 )=. 4b 
Pw  (4,4  )=—  1 . 3 
PW<  4, 5)= 1 .0 
PW (4,6) =0.0 
PN<5,  1  )=0.0 
Prt(b,2)=0.(. 

PW( 5, 3) =0.0 
Pri ( 5,4 )=0 .3 
Pri( 5,5)=-  1.  1  5 
PW( 5,6)=l .0 
PW(6, I )=0 .0 
PW C  6, 2 )=0.0 
PW( 6,3 )=0.0 
Prt(6,4)=0.0 
PW( 6,5 )=.  15 
PW(6,6)*-1.0 
GOTO  10 
5  CONTINUE 

DY(  I  )=-.6*Yp( I , I )+.5*YP(  I  ,2) 

0YC2.)=  ,6*YPi  I ,  I  >-l .  1*YP(  I  ,2)  +  YP(  i  ,3) 
DY(3)=.6*YP( I ,2 )“ I . 45*YP( 1 ,3>+YP( 1,4) 
DY(4)  =  .45*YP( 1 ,3)- 1 .3*YP( 1 ,4)+YP(  1,5) 
DY(5 )= .  3*YP( I ,4)-1 . 15*YP( 1 ,5)+YP( I ,6) 
DY(6)=. 15*YP( 1 ,5)-YP( 1,6) 

10  RETURN 
END 


Figure  A3.1 — continued. 
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To  cope  with  the  expanding  technology,  our  society  must 
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Applied  Science  is  completely  committed  to  this  ob¬ 
jective. 


